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ABSTRACT 

The  National  Bureau  of  Standards  (NBS),  Boulder,  Colorado,  uses  an  algorithm  which 
generates  UTC(NBS)  from  its  ensemble  of  clocks  and  automatically,  optimally,  and 
dynamically  weights  each  clock  in  the  ensemble.  The  same  algorithm  was  used  at  the 
Master  Control  Station  (MCS)  of  the  Global  Positioning  System  (GPS)  to  generate 
a  time  scale  from  a  small  ensemble  of  cesium  clocks.  Time  transfer  employing  the 
GPS  common  view  technique  between  NBS  (Boulder)  and  MCS  (Colorado  Springs) 
was  used  to  evaluate  the  stability  of  the  MCS  ensemble  relative  to  UTC(NBS).  The 
results  demonstrate  the  power  of  the  NBS  algorithm  in  providing  a  stable  time  scale 
from  a  small  ensemble  of  clocks.  The  resulting  scale  is,  in  principle,  more  stable  than 
the  best  clock  and  a  poor  clock  need  not  degrade  the  ensemble. 

I.  GENERAL 

We  decided  to  experimentally  incorporate  the  National  Bureau  of  Standards’  (NBS)  algorithm  to  generate 
ensemble  time  from  a  small  ensemble  of  commercial  cesium  beam  clocks  at  the  GPS  Master  Control  Station 
(MCS)  in  Colorado  Springs.  We  therefore  considered  using  GPS  common-view  time  transfer  between  the 
MCS  and  NBS,  Boulder  Colorado,  to  evaluate  the  result.  With  this  technique,  independently  generated  time 
scales  can  be  compared  over  long  distances  on  a  routine  daily  basis.  In  our  experiment,  the  distance  was 
only  90  miles,  so  differential  ephemeris  and  propagation  path  errors  were  practically  zero  in  the  common- 
view  mode  of  time  transfer.  Therefore,  daily  comparisons  of  MCS  ensemble  time  with  UTC(NBS)  were 
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accomplished  with  a  precision  of  1.4  nanoseconds,  providing  a  high  resolution  measure  of  the  MCS  ensemble 
performance. 

The  NBS  algorithm  generates  ensemble  time  from  a  collection  of  cesium  clocks;  in  place  in  the  Boulder 
laboratory  since  1968,  the  algorithm  produces  a  time  scale  that  is  the  weighted  mean  of  all  clocks  in  the 
ensemble.  Weight  assignment  is  semi-automatic  and  dynamic.  The  squared  error  of  the  scale  (defined  with 
respect  to  a  perfect  time  scale)  is  minimized. 

For  this  experiment,  the  NBS  time  scale  was  considered  to  be  perfect.  The  ensemble  at  the  MCS,  on  the  other 
hand,  had  lower  quality,  consisting  of  only  four  commercial  Hewlett-Packard  cesium  beam  clocks  (model 
5061 A),  one  of  which  was  not  a  high-performance  unit  (004  option).  Additionally,  no  special  environmental 
controls  were  in  place  other  than  the  normal  building  air  conditioning.  Nevertheless,  the  time  scale  produced 
by  the  MCS  ensemble  operating  with  the  NBS  algorithm  was  remarkably  stable  over  the  eighty  day  period 
of  the  experiment.  This  can  be  attributed  to  the  relatively  good  performance  of  two  clocks  in  the  ensemble 
and  the  fact  that  the  algorithm  automatically  and  appropriately  assigned  the  higher  weights  to  these  two 
clocks. 

II.  NBS  TIMING  ENSEMBLE  ALGORITHM 

The  following  is  an  outline  of  the  NBS  algorithm.  A  detailed  discussion  of  this  algorithm  can  be  found  in 
reference  2. 

As  input,  the  algorithm  requires  periodic  time-difference  measurements  between  all  pairs  of  clocks  in  the 
ensemble,  initial  values  of  time  and  frequency  for  each  clock,  and  two  parameters  derived  from  the  observed 
Allan  Variance  of  each  clock  (calculated  from  the  r-1/2  behavior  and  from  the  ”t”  of  the  Allan  Variance  at 
the  minimum  variance  point). 

As  output,  the  algorithm  produces  an  estimate  of  the  time  and  frequency  offset  of  each  clock  with  respect  to 
the  weighted  mean  of  all  clocks  (ensemble  time)  at  each  measurement  time.  Physical  realization  of  ensemble 
time  can  be  had  by  physically  steering  any  clock  in  the  ensemble  to  produce  an  average  zero-time  estimate  for 
the  steered  clock.  In  this  experiment,  no  steering  was  done,  but  one  of  the  ensemble  clocks  provided  timing 
for  the  GPS  time  transfer  receiver.  Through  this  clock  the  MCS  ensemble  was  compared  to  UTC(NBS). 

The  equations  and  definitions  are  as  follows: 

Definitions: 

Fi(t)  =  estimates  of  time  and  frequency  offsets  respectively  of  clock  i  at  time  t  with  respect  to  some 
reference  time  scale, 

Xi(t)  =  predicted  time  offset  of  clock  i  at  time  t, 

Yi(t)  =  estimate  of  frequency  of  clock  i  at  t  over  the  interval  t  to  t +t, 

Xult)  —  measured  time  difference  between  clocks  i  and  j  at  time  t, 

£i(g)  =  accumulated  error  in  time  estimate  of  clock  i  over  the  interval  r, 

(e|(r))  =  mean  squared  error  in  ensemble  time  over  the  interval  r  at  time  t, 

(  )  =  indicates  time  average, 
r  =  time  interval  between  measurements, 

Nr  —  time  constant  of  exponential  filter  to  estimate  the  current  mean  squared  error, 
n  =  number  of  clocks  in  the  ensemble, 

TMiNi  —  value  of  r  at  minimum  cy(T)  on  Allan  Variance  curve  for  clock  i. 
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Equations: 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 


(10) 
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Explanation  of  equations: 

Equation  1:  forms  a  prediction  of  the  time  offset  for  each  clock  for  the  next  measurement  time  (t+r)  based 
on  the  current  estimates  of  time  and  filtered  frequency. 

Equation  2:  estimates  the  time  offset  of  each  clock  j  at  time  t+r  given  the  measurements  X ij(t  +  r). 

Equation  3:  is  an  estimate  of  the  average  frequency  of  each  clock  over  the  interval  r  based  on  the  latest  two 
estimates  of  X,-. 

Equation  4:  incorporates  past  measurements  into  an  exponential  filtered  estimate  of  the  current  average 
frequency  of  clock  i.  The  time  constant  for  the  filter  was  chosen  to  be  NT  =  20  days.  The  exponential 
frequency-weighting  time  constant  (m,)  is  determined  from  the  relative  levels  of  white  noise  and  random 
walk  (or  flicker)  FM  for  clock  i  (equation  10). 


15 


Equation  5:  is  the  accumulated  error  in  the  estimate  of  X{  over  the  interval  r.  The  additive  term  K{  accounts 
for  the  fact  that  the  term  in  brackets  on  the  right-hand  side  of  Equation  5  is  biased  because  clock  i  is  part 
of  the  ensemble.  See  equation  9  to  calculate  K 

Equation  6:  is  an  exponential  time  filter  for  the  determination  of  the  mean  square  time  error  of  each  clock. 
Recognizing  that  the  noise  characteristics  of  a  cesium  clock  may  not  be  stationary,  past  measurements  are 
de-weighted  in  the  averaging  process.  The  initial  value  of  (e2(r))  can  be  estimated  as  r2er2(r). 

Equation  7:  forms  an  estimate  of  ensemble  time  error.  Any  clock  can  only  improve  this  number — a  poorly 
performing  clock  cannot  harm  the  stability  of  the  ensemble. 

Equation  8:  calculates  the  weight  to  be  used  in  Equation  2  for  each  clock.  When  calculated  this  way,  the 
resulting  error  in  ensemble  time  with  respect  to  a  perfect  clock  can  be  shown  to  be  minimized  in  a  least 
squares  sense. 

Equation  9:  The  error  estimate  from  the  first  term  on  the  right  of  Equation  5  is,  on  the  average,  biased 
small,  because  each  clock  is  a  member  of  the  ensemble,  and  sees  itself  through  its  weighting  factor.  The 
larger  a  clock’s  weight,  the  larger  is  the  bias.  Under  the  assumption  of  a  normal  distribution  of  errors  the 
size  of  the  bias  can  be  estimated  as  given  by  Equation  9,  which  is  added  to  Equation  5  in  order  to  remove 
the  bias,  on  the  average. 

Equation  10:  Computes  ny  used  in  Equation  4  to  form  the  filtered  estimate  of  the  frequency  of  clock  i.  This 
value  of  M  can  be  shown  to  minimize  the  error  in  predicting  time  (Equation  l)  given  two  kinds  of  noise  in 
the  clock  (white  and  random  walk  FM).  If  white  FM  and  flicker  FM  are  more  suitable  models,  then  mt-  can 
be  approximated  as  r//r,  where  r/  is  the  intercept  value  of  r  on  a  cry (r)  plot  for  the  white  and  flicker  FM. 

III.  DESCRIPTION  OE  EXPERIMENT 

Figure  1  is  a  combined  hardware/software  flow  chart  of  the  experiment.  Four  clocks,  located  at  the  MCS, 
are  linked  by  a  measurement  system  to  provide  the  time  difference  between  each  clock  and  a  reference  clock. 
Measurements  are  taken  hourly,  but,  for  this  experiment,  daily  measurements  were  used  to  eliminate  the 
effects  of  measurement  noise.  Differences  between  measurements  provide  the  data  required  for  the  NBS 
algorithm.  The  algorithm,  implemented  at  the  MCS  on  a  PC,  produced  the  ensemble  time  of  each  clock 
(ENS-1,  ENS-2,  etc.). 

The  reference  clock  also  provided  receiver  timing  for  common-view  time  transfer  between  MCS  and  NBS. 
The  common  view  time  transfer  data  were  processed  at  the  NBS  Time  and  Frequency  laboratory  to  provide 
daily  measures  of  the  time  offset  between  UTC(NBS)  and  the  reference  clock  in  the  MCS  ensemble  (NBS- 
REF).  These  data  were  differenced  with  the  measurements  (REF-1,  REF-2,  etc.)  so  the  time  behavior  of  each 
clock  in  the  MCS  ensemble  could  be  evaluated  with  respect  to  UTC(NBS).  Finally,  the  ensemble  behavior 
(i.e.  ENS-UTC(NBS))  was  obtained  from  each  clock’s  behavior  with  respect  to  both  the  MCS  ensemble  and 
UTC(NBS).  This  served  as  an  independent  measure  of  ensemble  behavior  over  the  test  period. 

IV.  RESULTS 

Figure  2  is  a  plot  of  the  accumulated  time  difference  between  MCS  ensemble  time  and  UTC  (NBS)  as  obtained 
in  the  manner  just  described.  A  constant  frequency  difference  of  1.2  parts  in  1013  has  been  removed.  A 
frequency  difference  will  exist  between  two  independent  time  scales;  the  fact  that  it  was  so  small  in  this  case 
is  coincidental.  The  nondeterministic  behavior  of  the  ensemble  is  quantified  in  the  Allan  Variance  curve 
(figure  3).  The  stability  of  the  four-clock  ensemble  is  remarkable  in  view  of  the  poor  performance  of  clocks  1 
and  2  (figure  5).  This  stability  is  attributed  to  the  good  performance  of  the  other  two  clocks,  and  the  NBS 
algorithm  which  weighted  the  better  clocks  higher  (figure  4). 

Figure  4  shows  the  history  of  the  individual  weights  assigned  by  the  algorithm  to  each  clock  over  the  eighty- 
day  period.  Clock  2  received  a  very  low  weight.  It  is  the  poorest  performing  clock  as  indicated  by  independent 
comparison  with  UTC(NBS)  (figure  5)  and  is  the  clock  without  the  high  performance  tube.  Because  of  the 
limited  time  span  of  the  test,  the  initial  weights  were  based  upon  an  experimental  run  with  the  first  two 
weeks’  data.  For  the  experimental  run,  the  clocks  were  initially  weighted  equally  (25%)  and  the  resulting 
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weights,  after  two  weeks,  were  used  to  initialize  the  eighty-day  run.  This  procedure  eliminated  a  two-week 
transient  period  which  would  have  biased  the  result. 

Note  that,  in  Figure  4,  the  weight  of  clock  1  dropped  to  zero  on  MJD  46971.  This  was  the  result  of 
an  automatic  de-weighting  of  this  clock  due  to  a  failure  of  the  algorithm  to  predict  this  clock’s  behavior 
between  measurements  to  within  three  sigma  of  the  estimated  clock  error.  Apparently  the  frequency  change 
exhibited  by  clock  1  (figure  5)  on  MJD  46971  was  detected  and  the  algorithm  de-weighted  clock  1,  preventing 
perturbation  of  the  ensemble.  In  figure  2,  removal  of  a  heavily  weighted  clock  did  not  perturb  ensemble  time 
or  frequency.  It  is  a  very  important  feature  of  the  NBS  algorithm  that  clocks  can  be  removed  (and  added) 
with  minimal  perturbation. 

The  weights  of  the  clocks  shown  in  figure  4  are  sometimes  called  the  short-term  stability  weights  since 
these  are  the  optimum  weights  for  the  interval  between  measurements.  The  ratio  of  the  long-term  stability 
performance  to  the  short-term  often  varies  significantly  between  clocks.  The  m,-  parameters  reflect  the 
long-term  stability  performance.  The  values  used  for  this  parameter  for  clocks  1,  2,  3,  and  4,  were  0.15,  1, 
18,  and  2  respectively.  Figure  6  shows  the  corresponding  frequency  stability,  <ry[r),  of  each  of  the  clocks  in 
the  MCS  ensemble. 

In  order  to  make  an  estimate  of  the  measurement  uncertainty  between  the  MCS  and  UTC(USNO  MC) 
(direct  common  view  time  transfer  between  the  MCS  and  USNO  is  carried  out  routinely)  NBS  measured 
UTC(USNO  MC)  minus  UTC(NBS)  (Figure  7),  then  differenced  that  with  the  data  in  figure  2.  This  gave 
a  different  path  estimate  of  UTC(USNO  MC)  minus  MCS  ensemble  than  the  direct  one.  The  difference 
between  the  direct  path  measurement  and  this  latter  path  is  plotted  in  figure  8,  and  the  corresponding  MOD 
CTj ,(r)  plot  is  shown  in  figure  9.  The  -1  slope  characterizing  the  stability  values  plotted  in  figure  9  would 
indicate  a  model  of  flicker  noise  time  (phase)  modulation  (PM)  with  Toy{r)  =  4  nanoseconds. 

V.  CONCLUSION 

The  theory  upon  which  the  NBS  timekeeping  algorithm  is  based  predicts  an  ensemble  stability  better  than 
the  best  clock  and,  in  principle,  the  ensemble  stability  is  not  harmed  by  a  poor  clock.  Also,  the  ensemble 
stability  is  not  adversely  perturbed  by  adding  or  dropping  a  clock.  This  experiment  supports  these  theoretical 
inferences;  the  experiment  is  based  on  real  data  and  independent  measures  of  ensemble  behavior. 

This  improved  clock  performance  from  a  weighted  clock  set,  when  employed  in  GPS  operationally,  will 
provide  three  significant  advantages:  (1)  the  improved  stability  will  provide  the  opportunity  to  have  better 
synchronization  to  UTC;  (2)  the  system  will  be  much  more  immune  to  the  failure  or  mal-performance  of 
an  individual  clock;  and  (3)  the  improved  long-term  performance  gained  from  the  ensemble  will  improve 
autonomous  GPS  performance  (i.e.  180  day  navigation  messages). 
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Figure  1.  A  combined  flow  chart  and  block  diagram  of  the  hardware/software  system 
employed  in  the  experiment. 
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Figure  2.  A  plot  of  the  residual  time  differences  between  the  MCS  ensemble  time  and 
UTC(NBS)  after  subtracting  the  mean  frequency  and  the  mean  time  from  the  data. 
The  plotted  values  cover  the  period  from  21  May  through  9  August  1987. 

Figure  3.  A  plot  of  the  fraction  frequency  stability,  uy{r),  of  the  data  shown  in  Figure 
2.  The  outstanding  long  term  stability  at  integration  time  of  a  few  weeks  is  evident. 
The  stability  at  one  and  two  days  is  probably  limited  by  measurement  noise. 
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Figure  4.  A  plot  of  the  percentage  weights  assigned  to  each  of  the  four  clocks  over  the 
course  of  the  experiment.  Notice  that  the  weight  of  clock  1  is  set  to  zero  at  the  point 
where  a  frequency  step  was  detected  (MJD  46970). 
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Figure  5.  A  plot  of  the  time  residuals  of  each  of  the  clocks  versus  the  MCS  ensemble 
after  removing  a  mean  frequency  and  also  after  setting  an  arbitrary  scaling  factor  for 
the  initial  ordinate  for  each  clock.  This  was  done  for  the  convenience  of  plotting  and 
display  of  these  residuals. 
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Figure  6.  A  plot  of  the  fractional  frequency  stability,  cry(r),  for  each  of  the  MCS  clocks 
against  an  independent  reference,  UTC(NBS).  The  advantage  of  an  algorithm  using 
statistical  weighting  is  obvious  from  the  very  different  performance  of  these  four  clocks. 
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Figure  7.  A  plot  of  UTC(USNO  MC)  versus  UTC(NBS)  is  shown  for  reference  for  the 
period  of  the  experiment.  Since  these  are  independent  time  scales,  we  can  conclude 
that  the  stability  shown  is  the  square  root  of  the  sum  of  the  variances  of  each,  and  that 
the  reference  scale  used  in  this  experiment,  UTC(NBS),  is  everywhere  better  than  the 
values  plotted  here. 
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Figure  8  &:  9.  A  test  was  made  on  the  GPS  common-view  measurement  residuals 
by  completing  the  measurement  triangle:  MCS  at  Colorado  Springs,  CO  to  NBS  at 
Boulder,  CO  to  USNO  at  Washington  D.  C.  and  thence  back  to  the  MCS.  Figure  8 
shows  the  time  residuals  for  the  sum  of  the  three  legs,  which  should  be  zero  if  the 
system  were  perfect,  and  Figure  9  is  a  modified  ay(r)  plot  showing  that  the  residuals 
can  be  modeled  by  a  flicker  noise  phase  modulation  (PM)  at  a  level  of  rav(r)  =  4 
nanoseconds  with  the  data  being  taken  once  per  sidereal  day  for  the  GPS  common- 
view  measurements.  If  each  leg  were  equal  in  its  contribution  and  independent  of  the 
others,  its  contribution  would  be  2.3  nanoseconds.  We  measured  the  common-view 
noise  of  the  MCS,  NBS  path  at  1.4  nanoseconds,  which  is  reasonable  since  that  leg  is 
much  shorter  than  the  other  two.  Notice  that  there  is  an  apparent  bias  of  about  -5 
nanoseconds  indicated  in  Figure  8. 
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QUESTIONS  AND  ANSWERS 

Gernot  Winkler,  USNO  (not  into  the  microphone):  Indecipherable  question 
Mr.  Varnum:  None  of  them.  That  clock  was  not  included  in  the  ensemble. 

Dr.  Winkler  had  a  further  question 

Mr.  Varnum:  Two  of  them  were  the  observatory,  two  of  them  were  GPS  clocks. 

Don  Percival,  University  of  Washington:  The  one  clock  which  experienced  the  abnormality  in  frequency 
and  seemed  to  recover  afterwards,  but  the  weight  did  not  start  to  increase  after  a  while.  Do  you  have  an 
explanation  for  that? 

Mr.  Varnum:  It  did  increase  but  it  was  not  visible  with  the  resolution  on  the  plot.  Perhaps  it  did  not 
increase  as  fast  as  it  should  have,  and  that  is  probably  a  limitation  within  the  algorithm  as  we  implemented 
it.  It  did  perform  well  enough  to  have  been  more  heavily  weighted  than  it  was. 

Mr.  Percival:  Did  you  give  any  thought  to  looking  at  the  Jones-Trion(sp?)  algorithm  which  was  developed 
at  NBS  to  see  how  it  would  compare  with  the  one  you  have  now? 

Mr.  Varnum:  No,  I  haven’t,  I  am  not  familiar  with  it. 

David  Allan,  NBS:  The  Jones-Trion{sp?)  that  Don  is  referring  to  is  based  on  Kalman  optimization.  In 
theory  it  is  supposed  to  be  as  good  or  better  than  the  NBS  algorithm.  We  have  found,  in  our  implementation 
of  it,  that  we  have  seen  frequency  steps  occur  when  clocks  are  added  or  deleted  from  the  ensemble.  We  are 
not  sure  why,  there  may  be  some  theoretical  problems.  The  NBS  algorithm  seems  to  be  immune  from  that 
effect  as  this  paper  demonstrates. 

Jim  Wright,  Pan  American:  How  long  did  it  take  you  to  establish  your  initial  parameters,  initial 
weighting? 

Mr.  Varnum:  I  used  a  two  week  period  for  that.  In  fact,  to  give  you  the  details,  I  processed  the  data  that 
you  saw  for  two  weeks  to  generate  a  set  of  weights,  then  went  back  to  the  beginning  and  used  those  weights 
as  the  input  to  the  algorithm. 

Mr.  Wright:  Another  question-when  you  dropped  down  to  three  clocks,  it  seems  that  one  of  the  clocks 
was  carrying  about  sixty  per  cent  of  the  entire  weight.  Do  you  have  any  concerns  about  that? 

Mr.  Varnum:  Yes,  that  is  a  concern  in  a  practical  case.  If  any  one  clock  assumes  too  much  weight  and 
something  happens  to  that  clock,  then  you  would  expect  that  the  ensemble  time  would  deteriorate.  In 
practice,  there  is  usually  a  limitation  applied.  I  did  not  do  that  in  this  experiment,  but  you  would  normally 
try  to  prevent  that  from  happening. 
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